#!/usr/bin/env python
# -*- coding: utf-8 -*-
import numpy as np


def eamgen(para):
    Nr = 10000
    rk = []
    ak = []
    
    # 确保para是有效数组
    if para is not None and hasattr(para, '__len__') and len(para) > 0:
        npairs = int(len(para) / 2)
        rk = para[:npairs]
        ak = para[npairs:]
    
    # 添加错误检查和默认值
    if len(rk) == 0 or len(ak) == 0 or len(rk) != len(ak):
        print("警告: 输入参数无效或长度不一致，使用默认参数")
        # 设置默认参数以确保生成有效的势能表
        rk = np.array([4.0])  # 默认截止距离
        ak = np.array([0.1])  # 默认振幅参数
    
    # 确保rk中只有正数
    rk = np.array([abs(r) for r in rk])
    
    out = []
    kw = "Re_He"
    cutoff = max(rk)
    
    # 确保cutoff是合理的正值
    if cutoff <= 0:
        cutoff = 4.0  # 默认截止距离
        print("警告: 计算的截断距离无效，使用默认值4.0")
    
    dr = cutoff / Nr  # 步长
    out.append(f"# Pair potential, Nr={Nr},  dr={dr:.12e}, cutoff={cutoff:.12e}\n\n{kw}\n")
    out.append(f"N  {Nr}  R   {dr:.12e}   {cutoff:.12e}\n\n")
    for i in range(1, Nr + 1):
        r = i * dr
        y = 0   # 势能
        dy = 0  # y一阶导数
        for j in range(len(rk)):
            y += ak[j] * (rk[j] - r) ** 3 * np.heaviside(rk[j] - r, 1.0)
            dy += -3 * ak[j] * (rk[j] - r) ** 2 * np.heaviside(rk[j] - r, 1.0)
        out.append(f"{i:<8d}   {r:.12e}  {y:.12e}  {-dy:.12e}\n")

    file = open("pot_Re_He", "w", encoding="utf-8")
    file.writelines(out)
    file.flush()
    file.close()
    return


if __name__ == "__main__":
    para = np.loadtxt("eam", dtype=float, skiprows=0, unpack=True, usecols=(1))
    eamgen(para)
